#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
#       Kolmogorov-Smirnov statistical test based on the equation 3 of Decin et al (2004) and eq 5 of Decin et al (2000)
#       
#       Copyright 2011 Jeffrey Simpson <jeffrey@isolaptop>
#       
#       This program is free software; you can redistribute it and/or modify
#       it under the terms of the GNU General Public License as published by
#       the Free Software Foundation; either version 2 of the License, or
#       (at your option) any later version.
#       
#       This program is distributed in the hope that it will be useful,
#       but WITHOUT ANY WARRANTY; without even the implied warranty of
#       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#       GNU General Public License for more details.
#       
#       You should have received a copy of the GNU General Public License
#       along with this program; if not, write to the Free Software
#       Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
#       MA 02110-1301, USA.
#       
#       

#       The first part of this finds the ratio between the observed spectrum and the synthetic spectrum and sums this over the range specified

import sys, glob, os, atpy, math, modelatm, param_finder_modules
import matplotlib.pyplot as plt
from numpy  import *

wavels='4500'
wavelf='4600'
resolution='2.3'


#This section loads the van Loon data and at the moment takes a given LEID and extracts the row
tbl_overall = atpy.Table('data/v07_WFI_2MASS.vot')
tbl_overall.sort('LEID')

params_found = open('params_v07_WFI_2MASS.out','w')

counter=0
for i in range(1):# range(len(tbl_overall)):
	using_BmV=0; no_B=0; flag="";
	BmV="";	Bmag=""; Vmag=""; M_V=""; VmK=""; temp_BmV=""; temp_VmK=""; log_g=""; temp=""; J_2MASS=""; K_2MASS=""
	best_parms=["","","","","",""]
	
	counter=counter+1; print counter
	
	LEID=tbl_overall.data['LEID'][i]
	RA=tbl_overall.data['RAJ2000'][i]
	Dec=tbl_overall.data['DEJ2000'][i]
	Vmag=tbl_overall.data['Vmag_B09'][i]
	Bmag=tbl_overall.data['Bmag_B09'][i]
	nS=tbl_overall.data['nS'][i]
	spFile=tbl_overall.data['spFile'][i]
	FeH_me=-1.49999
	params_found_KS_array=zeros( (192,7) )  
	params_found_chi_array=zeros( (192,7) ) 
	print LEID
	
	BmV=(Bmag-Vmag)-0.12
	Vmag=Vmag-0.36 #This is the A(V)		
	temp_BmV=param_finder_modules.Alonso_BmV(BmV,FeH_me)
	if temp_BmV=="":
		flag=flag+"b"
	J_2MASS=tbl_overall.data['j_m'][i]
	K_2MASS=tbl_overall.data['k_m'][i]
	#The A(V) is unapplied since the 2MASS photometry has not been corrected
	[JmK,VmK]=param_finder_modules.TCS_transforms(J_2MASS,K_2MASS,Vmag+0.36)
	temp_VmK=param_finder_modules.Alonso_VmK(VmK,FeH_me)
	temp_JmK=param_finder_modules.Alonso_JmK(JmK,FeH_me)
	if temp_VmK=="":
		temp=temp_BmV
		using_BmV=1
		flag=flag+"d"
	else:
		using_BmV=0
		temp=temp_VmK
	M_V=Vmag-13.7
	if temp!="" and temp<5750:
		temp_round=int(round(temp/250)*250)
		log_g=param_finder_modules.Alonso_logg(temp,FeH_me,M_V)
		if log_g=="":
			flag=flag+"e"
			continue
		grav_round=int(round(log_g/0.5)*0.5)
		[wave, intensity, specnum]=param_finder_modules.load_vanloon_spectrum(nS,spFile)
		#Loads the synthetic spectrum
		model_number=0
		#Finds the index of the first wavelength point in synth that is in obs
		[wave_synth,intensity_synth]=param_finder_modules.load_model_spectrum(4000,0.5,'5',resolution,wavels,wavelf,'0.0')
		[I,K_start,K_end]=param_finder_modules.indexfinder(wavels,wavelf,wave,wave_synth)
		abundance=0.0
		for FeH in [-0.5, -1.0, -1.5, -2.0]:
			print temp_round, grav_round, FeH
			metallicity=str(math.trunc(FeH*-10))
			[wave_synth,intensity_synth]=param_finder_modules.load_model_spectrum(temp_round,grav_round,metallicity,resolution,wavels,wavelf,'0.0')
			avg_cont_shift=param_finder_modules.avgcont(wave,intensity,wave_synth,intensity_synth)
			for abundance in [0.0, 0.5, 1.0]:
				model_number=model_number+1
				#[wave_synth,intensity_synth]=param_finder_modules.load_model_spectrum(temp_round,grav_round,metallicity,resolution,wavels,wavelf,abundance)
				for index, object in enumerate(intensity_synth):
				    intensity_synth[index] = object+avg_cont_shift
				#Calculates the KS beta parameter per Decin et al (2004) and (2000). 1<=k<=n-1 hence the ranges that are specified in the loops
				#I have tested this to make sure it is giving the right numbers.
				sup_trials=0.0
				sup_trials_best=0.0
				for k in range(len(I)-1):
					#This is the denominator of the ratio for eq 3
					V_k_denom=0.0
					for t in range(len(I)):
						V_k_denom=V_k_denom+intensity[K_start+t]/intensity_synth[I[t]]
					##This is the numerator of the ratio for eq 3
					V_k_numen=0.0
					for t in range(k+1):
						V_k_numen=V_k_numen+intensity[K_start+t]/intensity_synth[I[t]]
					sup_trials=abs(V_k_numen/V_k_denom-float(k+1)/float(len(I)))
					if sup_trials>sup_trials_best:
						sup_trials_best=sup_trials
				params_found_KS_array[model_number-1]=([model_number,temp_round,grav_round,FeH,abundance,sup_trials_best,0])
				chi2_1=param_finder_modules.chi2(I,K_start,K_end,intensity,intensity_synth)
				params_found_chi_array[model_number-1]=([model_number,temp_round,grav_round,FeH,abundance,chi2_1,0])
		
		#Sorts and ranks the two arrays so that the ranking can be used later. They are then resorted by model number.
		params_found_KS_array=params_found_KS_array[params_found_KS_array[:,5].argsort(),]
		m=0
		for k in range(len(params_found_KS_array)):
			m=m+1
			params_found_KS_array[k,6]=m
		params_found_chi_array=params_found_chi_array[params_found_chi_array[:,5].argsort(),]
		m=0
		for k in range(len(params_found_chi_array)):
			m=m+1
			params_found_chi_array[k,6]=m
			
		params_found_KS_array=params_found_KS_array[params_found_KS_array[:,0].argsort(),]
		params_found_chi_array=params_found_chi_array[params_found_chi_array[:,0].argsort(),]
		
		
		##This section joins them together and prints them out as one big array.
		params_found_joined=hstack((params_found_KS_array,params_found_chi_array[:,5:7]))
		
		i=-1
		distance_beta_chi=zeros( (192,1) )
		for row in params_found_joined:
			i=i+1
			distance_beta_chi[i]=(sqrt(power(row[6], 2)+power(row[8], 2)))
		#print params_found_joined
		#print distance_beta_chi
		params_found_joined=hstack((params_found_joined,distance_beta_chi))
		params_found_joined=params_found_joined[params_found_joined[:,9].argsort(),]
		
		params_found.write(str(LEID)+','+str(RA)+','+str(Dec)+','+str(Bmag)+','+str(Vmag)+','+str(BmV)+','+str(J_2MASS)+','+str(K_2MASS)+','+str(VmK)+','+str(temp_BmV)+','+str(temp_VmK)+','+str(log_g))
		#for i in range(9):
			#params_found.write(','+str(params_found_joined[0,i+1]))
		params_found.write(','+str(nS)+","+flag+'\n')
